void ViewTopoEvery16x16()

{


    //fileData = new TFile("MergedOutput.root");
    fileData = new TFile("./Pythia/MergedOutput.root");

//    fileMC = new TFile("./Pythia/MergedOutput.root");

    //TList *listData = fileData->Get("outputPhi256");

    TList *listData = fileData->Get("MCtrueOutputPhi256");

//    TList *listMC = fileMC->Get("outputPhi256");
//    TList *listMCtrue = fileMC->Get("MCtrueOutputPhi256");


    HistTopoQfQb = new TH2F("HistTopoQbQf", "<QfQb>", 16, -9*TMath::Pi()/16, 23*TMath::Pi()/16, 15, -0.05, 1.45); //PhiConf, EtaConf
    HistTopoQfQb->GetXaxis()->SetTitle("#phi sep");
    HistTopoQfQb->GetYaxis()->SetTitle("#eta sep");

    HistTopoQb = new TH2F("HistTopoQb", "<Qb>", 16, -9*TMath::Pi()/16, 23*TMath::Pi()/16, 15, -0.05, 1.45); //PhiConf, EtaConf
    HistTopoQb->GetXaxis()->SetTitle("#phi sep");
    HistTopoQb->GetYaxis()->SetTitle("#eta sep");

    HistTopoQf = new TH2F("HistTopoQf", "<Qf>", 16, -9*TMath::Pi()/16, 23*TMath::Pi()/16, 15, -0.05, 1.45); //PhiConf, EtaConf
    HistTopoQf->GetXaxis()->SetTitle("#phi sep");
    HistTopoQf->GetYaxis()->SetTitle("#eta sep");

    HistTopoEnum = new TH2F("HistTopoEnum", "<QbQf> -<Qb><Qf>", 16, -9*TMath::Pi()/16, 23*TMath::Pi()/16, 15, -0.05, 1.45); //PhiConf, EtaConf
    HistTopoEnum->GetXaxis()->SetTitle("#phi sep");
    HistTopoEnum->GetYaxis()->SetTitle("#eta sep");

    HistTopoQfQf = new TH2F("HistTopoQfQf", "<Qf><Qf>", 16, -9*TMath::Pi()/16, 23*TMath::Pi()/16, 15, -0.05, 1.45); //PhiConf, EtaConf
    HistTopoQfQf->GetXaxis()->SetTitle("#phi sep");
    HistTopoQfQf->GetYaxis()->SetTitle("#eta sep");

    HistTopoQf2 = new TH2F("HistTopoQf2", "<Qf*Qf>", 16, -9*TMath::Pi()/16, 23*TMath::Pi()/16, 15, -0.05, 1.45); //PhiConf, EtaConf
    HistTopoQf2->GetXaxis()->SetTitle("#phi sep");
    HistTopoQf2->GetYaxis()->SetTitle("#eta sep");

    HistTopoDenom = new TH2F("HistTopoDenom", "<QfQf> -<Qf><Qf>", 16, -9*TMath::Pi()/16, 23*TMath::Pi()/16, 15, -0.05, 1.45); //PhiConf, EtaConf
    HistTopoDenom->GetXaxis()->SetTitle("#phi sep");
    HistTopoDenom->GetYaxis()->SetTitle("#eta sep");

    HistTopoB = new TH2F("HistTopoB", "<QbQf> - <Qb><Qf> / <QfQf> - <Qf><Qf>", 16, -9*TMath::Pi()/16, 23*TMath::Pi()/16, 15, -0.05, 1.45); //PhiConf, EtaConf
    HistTopoB->GetXaxis()->SetTitle("#phi sep");
    HistTopoB->GetYaxis()->SetTitle("#eta sep");


    for (Int_t EtaConf=0; EtaConf<=14; EtaConf++)
        for (Int_t PhiConf=1; PhiConf<=12; PhiConf++)

        {
            THnSparseF* SparseData = (THnSparseF*) listData->FindObject(Form("fPhi1%dEta%d",PhiConf,EtaConf));

            TH1D *Qb = SparseData->Projection(6,"E");
            TH1D *Qf = SparseData->Projection(7,"E");
            TH1D *QfQb = SparseData->Projection(8,"E");
            TH1D *QfQf = SparseData->Projection(9,"E");


            Double_t QfMean = Qf->GetMean();
            Double_t QbMean = Qb->GetMean();
            Double_t QfQbMean = QfQb->GetMean();
            Double_t QfQfMean = QfQf->GetMean();


//            THnSparseF* SparseMC= (THnSparseF*) listMC->FindObject(Form("fPhi1%dEta%d",PhiConf,EtaConf));

//            TH1D *NpfMCESD = SparseMC->Projection(1,"E");
//            TH1D *NmfMCESD = SparseMC->Projection(3,"E");
//            TH1D *QfMCESD = SparseMC->Projection(7,"E");

//            Double_t NpfMeanMCESD = NpfMCESD->GetMean();
//            Double_t NmfMeanMCESD = NmfMCESD->GetMean();
//            Double_t QfMeanMCESD = QfMCESD->GetMean();


//            THnSparseF* SparseMCtrue = (THnSparseF*) listMCtrue->FindObject(Form("fPhi1%dEta%d",PhiConf,EtaConf));

//            TH1D *NpfMCtrue = SparseMCtrue->Projection(1,"E");
//            TH1D *NmfMCtrue = SparseMCtrue->Projection(3,"E");
//            TH1D *QfMCtrue = SparseMCtrue->Projection(7,"E");

//            Double_t NpfMeanMCtrue = NpfMCtrue->GetMean();
//            Double_t NmfMeanMCtrue = NmfMCtrue->GetMean();
//            Double_t QfMeanMCtrue = QfMCtrue->GetMean();

//            if(!QfMeanMCESD || !NpfMeanMCESD || !NmfMeanMCESD){cout<<"error! line~ 58";return;}

            HistTopoQfQb->SetBinContent(PhiConf+4, EtaConf+1, QfQbMean);
            HistTopoQb->SetBinContent(PhiConf+4, EtaConf+1, QbMean);
            HistTopoQf->SetBinContent(PhiConf+4, EtaConf+1, QfMean);
            HistTopoEnum->SetBinContent(PhiConf+4, EtaConf+1, QfQbMean - QfMean*QbMean);
            HistTopoQfQf->SetBinContent(PhiConf+4, EtaConf+1, QfMean*QfMean);
            HistTopoQf2->SetBinContent(PhiConf+4, EtaConf+1, QfQfMean);
            HistTopoDenom->SetBinContent(PhiConf+4, EtaConf+1, QfQfMean - QfMean*QfMean);
            HistTopoB->SetBinContent(PhiConf+4, EtaConf+1, (QfQbMean - QfMean*QbMean)/(QfQfMean - QfMean*QfMean));


        }



    for (Int_t EtaConf=0; EtaConf<=14; EtaConf++)
        for (Int_t PhiConf=13; PhiConf<=16; PhiConf++)

        {
            THnSparseF* SparseData = (THnSparseF*) listData->FindObject(Form("fPhi1%dEta%d",PhiConf,EtaConf));

            TH1D *Qb = SparseData->Projection(6,"E");
            TH1D *Qf = SparseData->Projection(7,"E");
            TH1D *QfQb = SparseData->Projection(8,"E");
            TH1D *QfQf = SparseData->Projection(9,"E");


            Double_t QfMean = Qf->GetMean();
            Double_t QbMean = Qb->GetMean();
            Double_t QfQbMean = QfQb->GetMean();
            Double_t QfQfMean = QfQf->GetMean();


//            THnSparseF* SparseMC= (THnSparseF*) listMC->FindObject(Form("fPhi1%dEta%d",PhiConf,EtaConf));

//            TH1D *NpfMCESD = SparseMC->Projection(1,"E");
//            TH1D *NmfMCESD = SparseMC->Projection(3,"E");
//            TH1D *QfMCESD = SparseMC->Projection(7,"E");

//            Double_t NpfMeanMCESD = NpfMCESD->GetMean();
//            Double_t NmfMeanMCESD = NmfMCESD->GetMean();
//            Double_t QfMeanMCESD = QfMCESD->GetMean();


//            THnSparseF* SparseMCtrue = (THnSparseF*) listMCtrue->FindObject(Form("fPhi1%dEta%d",PhiConf,EtaConf));

//            TH1D *NpfMCtrue = SparseMCtrue->Projection(1,"E");
//            TH1D *NmfMCtrue = SparseMCtrue->Projection(3,"E");
//            TH1D *QfMCtrue = SparseMCtrue->Projection(7,"E");

//            Double_t NpfMeanMCtrue = NpfMCtrue->GetMean();
//            Double_t NmfMeanMCtrue = NmfMCtrue->GetMean();
//            Double_t QfMeanMCtrue = QfMCtrue->GetMean();

//            if(!QfMeanMCESD || !NpfMeanMCESD || !NmfMeanMCESD){cout<<"error! line~ 58";return;}

            HistTopoQfQb->SetBinContent(PhiConf-12, EtaConf+1, QfQbMean);
            HistTopoQb->SetBinContent(PhiConf-12, EtaConf+1, QbMean);
            HistTopoQf->SetBinContent(PhiConf-12, EtaConf+1, QfMean);
            HistTopoEnum->SetBinContent(PhiConf-12, EtaConf+1, QfQbMean - QfMean*QbMean);
            HistTopoQfQf->SetBinContent(PhiConf-12, EtaConf+1, QfMean*QfMean);
            HistTopoQf2->SetBinContent(PhiConf-12, EtaConf+1, QfQfMean);
            HistTopoDenom->SetBinContent(PhiConf-12, EtaConf+1, QfQfMean - QfMean*QfMean);
            HistTopoB->SetBinContent(PhiConf-12, EtaConf+1, (QfQbMean - QfMean*QbMean)/(QfQfMean - QfMean*QfMean));


        }

    HistTopoQfQb->SetBinContent(4,1,HistTopoQfQb->GetBinContent(5,3));
    HistTopoQfQb->SetBinContent(4,2,HistTopoQfQb->GetBinContent(5,3));
    HistTopoQfQb->SetBinContent(5,1,HistTopoQfQb->GetBinContent(5,3));
    HistTopoQfQb->SetBinContent(5,2,HistTopoQfQb->GetBinContent(5,3));
    HistTopoQfQb->SetBinContent(6,1,HistTopoQfQb->GetBinContent(5,3));
    HistTopoQfQb->SetBinContent(6,2,HistTopoQfQb->GetBinContent(5,3));


    HistTopoEnum->SetBinContent(4,1,HistTopoEnum->GetBinContent(5,3));
    HistTopoEnum->SetBinContent(4,2,HistTopoEnum->GetBinContent(5,3));
    HistTopoEnum->SetBinContent(5,1,HistTopoEnum->GetBinContent(5,3));
    HistTopoEnum->SetBinContent(5,2,HistTopoEnum->GetBinContent(5,3));
    HistTopoEnum->SetBinContent(6,1,HistTopoEnum->GetBinContent(5,3));
    HistTopoEnum->SetBinContent(6,2,HistTopoEnum->GetBinContent(5,3));

    HistTopoB->SetBinContent(4,1,HistTopoB->GetBinContent(5,3));
    HistTopoB->SetBinContent(4,2,HistTopoB->GetBinContent(5,3));
    HistTopoB->SetBinContent(5,1,HistTopoB->GetBinContent(5,3));
    HistTopoB->SetBinContent(5,2,HistTopoB->GetBinContent(5,3));
    HistTopoB->SetBinContent(6,1,HistTopoB->GetBinContent(5,3));
    HistTopoB->SetBinContent(6,2,HistTopoB->GetBinContent(5,3));


    //inverting HistTopo:
//    for (Int_t EtaConf=1; EtaConf<=8; EtaConf++)
//        for (Int_t PhiConf=1; PhiConf<=8; PhiConf++)

//        {

//            HistTopo->SetBinContent(PhiConf, EtaConf, -1*(HistTopo->GetBinContent(PhiConf,EtaConf)));

//        }

    gROOT->SetStyle("Plain");
    gStyle->SetPalette(1);

    TCanvas *c1 = new TCanvas("c1","name",50,50,800,500);
    c1->Divide(4,2,0.001,0.001);
    c1->cd(1);
    HistTopoQf->GetZaxis()->SetRangeUser(0,0.001);
    HistTopoQf->Draw("SURF1");
    c1->cd(2);
    HistTopoQb->GetZaxis()->SetRangeUser(0,0.001);
    HistTopoQb->Draw("SURF1");
    c1->cd(3);
    HistTopoQfQb->GetZaxis()->SetRangeUser(-0.001,0);
    HistTopoQfQb->Draw("SURF1");
    c1->cd(4);
    HistTopoEnum->GetZaxis()->SetRangeUser(-0.001,0);
    HistTopoEnum->Draw("SURF1");
    c1->cd(5);
    HistTopoQfQf->GetZaxis()->SetRangeUser(0,0.001);
    HistTopoQfQf->Draw("SURF1");
    c1->cd(6);
    HistTopoQf2->GetZaxis()->SetRangeUser(0,0.01);
    HistTopoQf2->Draw("SURF1");
    c1->cd(7);
    HistTopoDenom->GetZaxis()->SetRangeUser(0,0.01);
    HistTopoDenom->Draw("SURF1");
    c1->cd(8);
    HistTopoB->GetZaxis()->SetRangeUser(-0.01,0);
    HistTopoB->Draw("SURF1");

    return;

}
